Gamma-Poisson Model

Regression

Count Data

GLM

An extension of the Poisson model for count data that exhibits overdispersion.

General Principles

To model the relationship between a count outcome variable and one or more independent variables with overdispersion 🛈, we can use the Negative Binomial model.

Considerations

Caution

We have the same considerations as for the Poisson model.
Overdispersion is handled because the Gamma-Poisson model assumes that each Poisson count observation has its own rate. This is an additional parameter specified in the model (in the code, it is log_days).

Example

Below is an example code snippet demonstrating a Bayesian Gamma-Poisson model using the Bayesian Inference (BI) package.

Code

from BI import bi, jnp

# Setup device ------------------------------------------------
m = bi(platform='cpu') # Import

# Import Data & Data Manipulation ------------------------------------------------
# Import
from importlib.resources import files
data_path =  m.load.sim_gamma_poisson(only_path=True)
m.data(data_path, sep=',') 


# Define model ------------------------------------------------
def model(log_days, monastery, y):
    a = m.dist.normal(0, 1, name = 'a', shape=(1,))
    b = m.dist.normal(0, 1, name = 'b', shape=(1,))
    phi = m.dist.exponential(1, name = 'phi', shape=(1,))
    mu = jnp.exp(log_days + a + b * monastery)
    Lambda =  m.dist.gamma(rate = mu*phi, concentration = phi, name = 'Lambda')
    m.dist.poisson(rate = Lambda, obs=y)
# Run MCMC ------------------------------------------------
m.fit(model) # Optimize model parameters through MCMC sampling

# Summary ------------------------------------------------
m.summary() # Get posterior distributions

jax.local_device_count 16

  0%|          | 0/1000 [00:00<?, ?it/s]warmup:   0%|          | 1/1000 [00:01<23:14,  1.40s/it, 1 steps of size 2.34e+00. acc. prob=0.00]warmup:   1%|          | 6/1000 [00:01<03:09,  5.24it/s, 511 steps of size 1.14e-02. acc. prob=0.43]warmup:   1%|          | 11/1000 [00:01<01:41,  9.79it/s, 511 steps of size 5.90e-03. acc. prob=0.59]warmup:   2%|▏         | 16/1000 [00:01<01:15, 13.05it/s, 1023 steps of size 4.14e-03. acc. prob=0.64]warmup:   2%|▏         | 19/1000 [00:02<01:06, 14.75it/s, 127 steps of size 2.52e-02. acc. prob=0.70] warmup:   2%|▏         | 22/1000 [00:02<01:05, 14.82it/s, 511 steps of size 8.18e-03. acc. prob=0.69]warmup:   2%|▎         | 25/1000 [00:02<01:03, 15.44it/s, 1023 steps of size 3.45e-03. acc. prob=0.69]warmup:   3%|▎         | 28/1000 [00:02<00:55, 17.38it/s, 127 steps of size 1.13e-02. acc. prob=0.71] warmup:   3%|▎         | 33/1000 [00:02<00:50, 19.00it/s, 1023 steps of size 3.34e-03. acc. prob=0.71]warmup:   4%|▎         | 36/1000 [00:02<00:47, 20.13it/s, 127 steps of size 6.20e-03. acc. prob=0.72] warmup:   4%|▍         | 40/1000 [00:02<00:41, 22.95it/s, 511 steps of size 4.09e-03. acc. prob=0.72]warmup:   4%|▍         | 44/1000 [00:03<00:36, 26.00it/s, 255 steps of size 6.81e-03. acc. prob=0.73]warmup:   5%|▍         | 48/1000 [00:03<00:34, 27.55it/s, 511 steps of size 3.45e-03. acc. prob=0.73]warmup:   5%|▌         | 51/1000 [00:03<00:35, 26.85it/s, 127 steps of size 7.90e-03. acc. prob=0.74]warmup:   6%|▌         | 57/1000 [00:03<00:28, 32.99it/s, 255 steps of size 1.09e-02. acc. prob=0.75]warmup:   6%|▌         | 61/1000 [00:03<00:29, 31.30it/s, 255 steps of size 9.98e-03. acc. prob=0.75]warmup:   7%|▋         | 66/1000 [00:03<00:26, 35.72it/s, 3 steps of size 2.28e-03. acc. prob=0.74]  warmup:   7%|▋         | 70/1000 [00:03<00:28, 32.25it/s, 127 steps of size 7.37e-03. acc. prob=0.75]warmup:   7%|▋         | 74/1000 [00:03<00:28, 32.47it/s, 255 steps of size 4.66e-03. acc. prob=0.75]warmup:   8%|▊         | 78/1000 [00:04<00:28, 32.56it/s, 511 steps of size 3.81e-03. acc. prob=0.75]warmup:   8%|▊         | 83/1000 [00:04<00:26, 34.92it/s, 255 steps of size 5.87e-03. acc. prob=0.75]warmup:   9%|▉         | 88/1000 [00:04<00:25, 36.24it/s, 255 steps of size 5.61e-03. acc. prob=0.76]warmup:   9%|▉         | 92/1000 [00:04<00:28, 31.40it/s, 255 steps of size 8.23e-03. acc. prob=0.76]warmup:  10%|▉         | 96/1000 [00:04<00:28, 31.73it/s, 255 steps of size 5.14e-03. acc. prob=0.76]warmup:  10%|█         | 100/1000 [00:04<00:27, 33.21it/s, 127 steps of size 1.01e-02. acc. prob=0.76]warmup:  11%|█         | 107/1000 [00:04<00:21, 41.37it/s, 127 steps of size 7.50e-02. acc. prob=0.76]warmup:  12%|█▏        | 121/1000 [00:04<00:13, 66.17it/s, 31 steps of size 2.47e-01. acc. prob=0.77] warmup:  14%|█▎        | 137/1000 [00:05<00:09, 90.97it/s, 63 steps of size 9.63e-02. acc. prob=0.77]warmup:  15%|█▌        | 152/1000 [00:05<00:07, 106.83it/s, 1 steps of size 3.70e-01. acc. prob=0.77]warmup:  17%|█▋        | 169/1000 [00:05<00:06, 124.47it/s, 63 steps of size 1.21e-01. acc. prob=0.77]warmup:  18%|█▊        | 185/1000 [00:05<00:06, 133.70it/s, 63 steps of size 1.41e-01. acc. prob=0.78]warmup:  20%|██        | 201/1000 [00:05<00:05, 140.32it/s, 31 steps of size 1.42e-01. acc. prob=0.78]warmup:  22%|██▏       | 218/1000 [00:05<00:05, 147.72it/s, 31 steps of size 1.64e-01. acc. prob=0.78]warmup:  24%|██▎       | 235/1000 [00:05<00:04, 153.60it/s, 63 steps of size 1.38e-01. acc. prob=0.78]warmup:  26%|██▌       | 255/1000 [00:05<00:04, 162.53it/s, 63 steps of size 1.29e-01. acc. prob=0.78]warmup:  27%|██▋       | 272/1000 [00:05<00:04, 155.39it/s, 15 steps of size 3.96e-01. acc. prob=0.78]warmup:  29%|██▉       | 288/1000 [00:06<00:04, 152.52it/s, 31 steps of size 1.53e-01. acc. prob=0.78]warmup:  31%|███       | 309/1000 [00:06<00:04, 166.23it/s, 31 steps of size 1.15e-01. acc. prob=0.78]warmup:  33%|███▎      | 326/1000 [00:06<00:04, 166.73it/s, 31 steps of size 1.77e-01. acc. prob=0.78]warmup:  34%|███▍      | 343/1000 [00:06<00:03, 165.18it/s, 31 steps of size 8.23e-02. acc. prob=0.78]warmup:  36%|███▌      | 360/1000 [00:06<00:04, 159.33it/s, 31 steps of size 1.34e-01. acc. prob=0.78]warmup:  38%|███▊      | 377/1000 [00:06<00:03, 162.19it/s, 31 steps of size 1.85e-01. acc. prob=0.78]warmup:  39%|███▉      | 394/1000 [00:06<00:03, 163.78it/s, 31 steps of size 1.88e-01. acc. prob=0.78]warmup:  41%|████      | 412/1000 [00:06<00:03, 167.00it/s, 31 steps of size 1.71e-01. acc. prob=0.78]warmup:  43%|████▎     | 429/1000 [00:06<00:03, 163.72it/s, 31 steps of size 8.71e-02. acc. prob=0.78]warmup:  45%|████▍     | 446/1000 [00:06<00:03, 157.95it/s, 31 steps of size 2.10e-01. acc. prob=0.79]warmup:  46%|████▌     | 462/1000 [00:07<00:03, 141.83it/s, 63 steps of size 1.54e-01. acc. prob=0.78]warmup:  48%|████▊     | 477/1000 [00:07<00:04, 130.61it/s, 127 steps of size 9.02e-02. acc. prob=0.78]warmup:  49%|████▉     | 491/1000 [00:07<00:04, 125.05it/s, 31 steps of size 1.67e-01. acc. prob=0.78] sample:  51%|█████     | 506/1000 [00:07<00:03, 130.60it/s, 31 steps of size 1.33e-01. acc. prob=0.90]sample:  52%|█████▏    | 523/1000 [00:07<00:03, 138.78it/s, 31 steps of size 1.33e-01. acc. prob=0.89]sample:  54%|█████▍    | 539/1000 [00:07<00:03, 143.16it/s, 31 steps of size 1.33e-01. acc. prob=0.89]sample:  56%|█████▌    | 555/1000 [00:07<00:03, 146.90it/s, 31 steps of size 1.33e-01. acc. prob=0.87]sample:  57%|█████▋    | 571/1000 [00:07<00:02, 149.84it/s, 31 steps of size 1.33e-01. acc. prob=0.87]sample:  59%|█████▊    | 587/1000 [00:08<00:02, 148.14it/s, 31 steps of size 1.33e-01. acc. prob=0.87]sample:  60%|██████    | 602/1000 [00:08<00:02, 148.45it/s, 31 steps of size 1.33e-01. acc. prob=0.87]sample:  62%|██████▏   | 617/1000 [00:08<00:02, 147.17it/s, 31 steps of size 1.33e-01. acc. prob=0.88]sample:  63%|██████▎   | 633/1000 [00:08<00:02, 148.14it/s, 31 steps of size 1.33e-01. acc. prob=0.88]sample:  65%|██████▍   | 649/1000 [00:08<00:02, 149.88it/s, 31 steps of size 1.33e-01. acc. prob=0.88]sample:  66%|██████▋   | 665/1000 [00:08<00:02, 150.29it/s, 31 steps of size 1.33e-01. acc. prob=0.89]sample:  68%|██████▊   | 681/1000 [00:08<00:02, 151.83it/s, 31 steps of size 1.33e-01. acc. prob=0.89]sample:  70%|██████▉   | 697/1000 [00:08<00:01, 152.75it/s, 31 steps of size 1.33e-01. acc. prob=0.89]sample:  71%|███████▏  | 713/1000 [00:08<00:01, 154.13it/s, 31 steps of size 1.33e-01. acc. prob=0.89]sample:  73%|███████▎  | 729/1000 [00:08<00:01, 153.94it/s, 31 steps of size 1.33e-01. acc. prob=0.89]sample:  74%|███████▍  | 745/1000 [00:09<00:01, 153.57it/s, 31 steps of size 1.33e-01. acc. prob=0.89]sample:  76%|███████▌  | 761/1000 [00:09<00:01, 154.41it/s, 31 steps of size 1.33e-01. acc. prob=0.89]sample:  78%|███████▊  | 777/1000 [00:09<00:01, 151.51it/s, 31 steps of size 1.33e-01. acc. prob=0.89]sample:  79%|███████▉  | 793/1000 [00:09<00:01, 152.97it/s, 31 steps of size 1.33e-01. acc. prob=0.89]sample:  81%|████████  | 809/1000 [00:09<00:01, 152.69it/s, 31 steps of size 1.33e-01. acc. prob=0.90]sample:  82%|████████▎ | 825/1000 [00:09<00:01, 152.66it/s, 31 steps of size 1.33e-01. acc. prob=0.90]sample:  84%|████████▍ | 841/1000 [00:09<00:01, 152.96it/s, 31 steps of size 1.33e-01. acc. prob=0.90]sample:  86%|████████▌ | 857/1000 [00:09<00:00, 147.99it/s, 31 steps of size 1.33e-01. acc. prob=0.90]sample:  87%|████████▋ | 873/1000 [00:09<00:00, 149.93it/s, 31 steps of size 1.33e-01. acc. prob=0.89]sample:  89%|████████▉ | 889/1000 [00:10<00:00, 150.06it/s, 31 steps of size 1.33e-01. acc. prob=0.89]sample:  90%|█████████ | 905/1000 [00:10<00:00, 150.95it/s, 31 steps of size 1.33e-01. acc. prob=0.89]sample:  92%|█████████▏| 921/1000 [00:10<00:00, 150.37it/s, 31 steps of size 1.33e-01. acc. prob=0.89]sample:  94%|█████████▎| 937/1000 [00:10<00:00, 148.06it/s, 31 steps of size 1.33e-01. acc. prob=0.89]sample:  95%|█████████▌| 952/1000 [00:10<00:00, 145.21it/s, 31 steps of size 1.33e-01. acc. prob=0.89]sample:  97%|█████████▋| 967/1000 [00:10<00:00, 146.03it/s, 31 steps of size 1.33e-01. acc. prob=0.89]sample:  98%|█████████▊| 983/1000 [00:10<00:00, 149.36it/s, 31 steps of size 1.33e-01. acc. prob=0.89]sample: 100%|█████████▉| 998/1000 [00:12<00:00, 24.45it/s, 31 steps of size 1.33e-01. acc. prob=0.89] sample: 100%|██████████| 1000/1000 [00:12<00:00, 79.79it/s, 31 steps of size 1.33e-01. acc. prob=0.89]
arviz - WARNING - Shape validation failed: input_shape: (1, 500), minimum_shape: (chains=2, draws=4)

	mean	sd	hdi_5.5%	hdi_94.5%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
Lambda[0]	1.39	0.36	0.83	1.91	0.01	0.02	537.43	299.02	NaN
Lambda[1]	1.76	0.40	1.13	2.34	0.01	0.02	713.23	285.08	NaN
Lambda[2]	1.57	0.42	0.92	2.17	0.01	0.02	826.04	397.99	NaN
Lambda[3]	1.47	0.36	0.86	1.97	0.01	0.02	850.98	414.15	NaN
Lambda[4]	1.53	0.38	0.95	2.14	0.01	0.02	630.33	428.36	NaN
...	...	...	...	...	...	...	...	...	...
Lambda[3398]	3.04	0.77	1.93	4.19	0.03	0.04	854.46	465.41	NaN
Lambda[3399]	2.86	0.73	1.76	3.95	0.03	0.04	341.11	200.88	NaN
a[0]	-0.42	0.01	-0.44	-0.40	0.00	0.00	51.13	111.55	NaN
b[0]	-2.78	0.03	-2.83	-2.73	0.00	0.00	63.60	188.73	NaN
phi[0]	16.10	1.92	13.30	19.15	0.55	0.14	12.83	101.44	NaN

3403 rows × 9 columns

library(BayesianInference)

# Setup platform------------------------------------------------
m=importBI(platform='cpu')

# Import data ------------------------------------------------
m$data(m$load$sim_gamma_poisson(only_path=T), sep = '')
m$data_to_model(list('log_days', 'monastery', 'y' )) # Send to model (convert to jax array)

# Define model ------------------------------------------------
model <- function(log_days, monastery, y){
  # Parameter prior distributions
  alpha = bi.dist.normal(0, 1, name='alpha', shape=c(1))
  beta = bi.dist.normal(0, 1, name='beta', shape=c(1))
  phi = bi.dist.exponential(1, name='phi', shape=c(1))
  mu = jnp$exp(log_days + alpha + beta * monastery)
  Lambda =  bi.dist.gamma(rate = mu*phi, concentration = phi, name = 'Lambda')
  # Likelihood
  bi.dist.poisson(rate=Lambda, obs=y)
}
# Run MCMC ------------------------------------------------
m$fit(model) # Optimize model parameters through MCMC sampling

# Summary ------------------------------------------------
m$summary() # Get posterior distributions

using BayesianInference

# Setup device------------------------------------------------
m = importBI(platform="cpu")

# Import Data & Data Manipulation ------------------------------------------------
# Import
data_path = m.load.sim_gamma_poisson(only_path = true)
m.data(data_path, sep=',')

# Define model ------------------------------------------------
@BI function model(log_days, monastery, y)
    a = m.dist.normal(0, 1, name = "a", shape=(1,))
    b = m.dist.normal(0, 1, name = "b", shape=(1,))
    phi = m.dist.exponential(1, name = "phi", shape=(1,))
    mu = jnp.exp(log_days + a + b * monastery)
    Lambda =  m.dist.gamma(rate = mu*phi, concentration = phi, name = "Lambda")
    m.dist.poisson(rate = Lambda, obs=y)
end

# Run mcmc ------------------------------------------------
m.fit(model)  # Optimize model parameters through MCMC sampling

# Summary ------------------------------------------------
m.summary() # Get posterior distributions

Mathematical Details

Bayesian model

In the Bayesian formulation, we define each parameter with priors 🛈. We can express the Bayesian regression model accounting for prior distributions as follows:

Y_i \sim \text{Poisson}(\lambda_i)

\lambda_i \sim \text{Gamma}(\mu_i \phi, \phi)

\log(\mu_i) = \text{rates}_i + \alpha + \beta X_i

\alpha \sim \text{Normal}(0,1)

\beta \sim \text{Normal}(0,1)

\phi \sim \text{Exponential}(1)

Where:

Y_i is the dependent variable for observation i.
\lambda_i is the rate parameter of the Poisson distribution for observation i, assuming that each Poisson count observation has its own rate_i.
\mu_i is the mean rate parameter.
\phi controls the level of overdispersion in the rates.
\alpha is the intercept term.
\beta is the regression coefficient.
X_i is the value of the predictor variable for observation i.

Notes

Note

We can apply multiple variables similarly as in chapter 2.
We can apply interaction terms similarly as in chapter 3.
We can apply categorical variables similarly as in chapter 4.